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The second order response functions and susceptibilities of finite temperature Bose-Einstein Con- 
densates (BEC) in a one dimensional harmonic trap driven by an external field that couples to the 
particle density are calculated by solving the time-dependent Hartree-Fock-Bogoliubov (TDHFB) 
equations. These provide additional insight into BEC dynamics, beyond the linear response regime. 
The results also apply to electron liquids in the Fractional Quantum Hall Effect (FQHE) regime 
which can be mapped onto an effective boson system coupled to a Chern-Simons gauge field. 

PACS numbers: 03.75.Fi 

I. INTRODUCTION 

We have recently studied an externally driven, finite temperature Bose-Einstein condensates (BEC) described by the 
time-dependent Hartree-Fock-Bogoliubov (TDHFB) equations X]- A systematic procedure was outlined for solving 
these equations perturbatively in the applied external field, and position-dependent linear response functions and 
susceptibilities were calculated. 

In this paper we extend this formalism to calculate the second order response function for the condensate, non- 
condensate density and non-condensate correlation. The linear response provides an adequate description of the 
system only for weak external perturbations. Otherwise, nonlinear effects contained in the higher order terms in the 
perturbation series may not be ignored. There are strong similarities between nonlinear optics and the dynamics 
of BEC owing to the interatomic interactions in Bose condensates. Previous work done in this area includes the 
demonstration of four wave mixing in zero temperature BEC using the Gross-Pitaevskii Equation (GPE) [2j]. In the 
present article we study nonlinear properties of BEC using the TDHFB framework. 

Numerous dynamical theories exist for finite temperature BEC that takes into account higher order collision pro- 
cesses, such as the time-dependent Bogoliubov-de Gennes equations H, the Hartree-Fock-Bogoliubov (HFB) the- 
ory |1 H IE 0. Quantum Kinetic Theory [1 II El EJ El El , and Stochastic methods El 111 EI El El The 
TDHFB theory is a self-consistent theory of BEC in the collisionless regime that progresses logically from the Gross- 
Pitaevskii Equation by taking into account higher order correlations of noncondensate operators. Although TDHFB 
neglects higher order correlations included in the various quantum kinetic theories, the TDHFB equations are valid 
at very low temperatures near zero, even down to the zero temperature limit, and are far simpler than the kinetic 
equations which can only be solved using approximations such as ZNG. Another attractive feature of TDHFB from 
a purely pragmatic point of view is that the Fermionic version of the theory has already been well-developed in 
Nuclear Physics [4| . We therefore work at the TDHFB level in this paper and our approach draws upon the analogy 
with the time-dependent Hartree-Fock (TDHF) formalism developed for nonlinear optical response of many electron 
systems 19]. 

In Section II we introduce the second order time and frequency domain response functions for an externally driven 
BEC. Numerical results are discussed in Sections III and IV for a condensate of 2000 atoms in a one dimensional 
harmonic trap. In Section V we show how this formalism may be applied for computing the second order response for 
j_j i a Fractional Quantum Hall system. In Section VI we conclude. Details of the derivation of the second order response 
C$ ' functions are given in the Appendix. 



II. NONLINEAR RESPONSE FUNCTION OF EXTERNALLY DRIVEN BEC 



We adopt the notation of Ref. [l| throughout the paper. The Hamiltonian describing the system of an externally 
driven, trapped atomic BEC is given by: 



H = H +H'(t), 



(1) 
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where 

Ha=Y2 "u">"j + ^Y1 V ijk ia\a\a k a h (2) 

ijkl 



and 



H'{t)=fiY t E ii {t)a\a j (3) 

ij 

with 

Ea{t)= [d?r4>*{r)V f {r,t)h{r)a\a r (4) 



The boson operators a| ( a* ) create ( annihilate ) a particle from a basis state with wave functions <f>i (r) . The single 
particle Hamiltonian H sp in Eq. J2J is diagonal if the basis state <f>i(r) is chosen to be the eigenstates of the trap, 
while the symmetrized interaction matrix elements, 

V ijkl = J d 3 rd 3 r'<f>*(r)^(r')V(r - r')0 fe (r')^(r), (5) 

describe the collision between the atoms, with V(r — r') being a general interatomic potential. H'it) describes the 
effect of a general external forc e Vj (r, t) on the condensate that mimics the mechanical force applied experimentally 
such as shaking of the trap [2JJ, |3(| . 

The dynamics of the system is calculated by solving the time-dependent Hartree-Fock-Bogoliubov (TDHFB) equa- 
tions for the condensate mean field, Zi — (a^), the non- condensate density pij — {d\a,j) — (al)(dj), and the non- 
condensate correlations Kij — — {a,i)(a,j). These are presented in Appendix ^ These nonlinear coupled 
equations are solved by an order by order expansion of the variables Zj, and Kij] at each order, the resulting 
equations to be solved become linear pj . It is found that the sequence of linear equations to be solved has the general 
form: 

ih d$ri(fi = £(n) ^ (n) ^ + A(n) (<) (g) 
at 

where we have denoted the set of n'th order variables zf" , p\" , and n[j as a 2N(2N + 1) x 1 column vector in the 

Liouville space notation where N is the number of basis states used, t/j^ = [z^ n \ z^ n '*, p^ n \ K^ n \ p^ n ^* , k,^*] t , and 
the 2N(2N+1) x 2iV(2iV+l) matrix C [n) is the n-th order Liouville operator obtained from the TDHFB equations Q. 
The matrices for all orders n > are identical i.e. = = ■ ■ ■ £( n > so that only the matrices £ < - - ) and 
are required to be calculated. £(°) and £W are presented in Appendix 151 The formal solution to Eq. © is: 



n 



f«)(t) = 1 f exp 
in J 

In the frequency domain, the solution to Eq. © takes the form 



\ in \t')dt'. (7) 



^ (n) H = ^7w A(n)(a;) ' (8) 

where ip^ n \uj) and X^(uj) are the Fourier transforms of ip^ n \t) and \( n \t) respectively. 
For the zero'th order (n = 0), we obtain the time- independent HFB equations (TIHFB) 

£(0)^(0)(t) = 0, (9) 

while for the first order (n = 1), the equation solved is Eq. Jfjjl with A^(i) being a 2N(2N + 1) x 1 vector 
calculated in Ref. 0; ((t) is also presented in Appendix iBl 

Once the n'th order solution to TDHFB is found, we can proceed to define the n'th order response functions. The 
physical significance of the response functions become more transparent when the n'th order solutions where a 
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is one of the variables z, p or k are expressed in real space. We therefore introduce the position dependent variables 
written in terms of the trap eigenstate basis: 



*W(r,t) = X>j B m(r), (10) 

3 

pW(r,t) = ^p^(*)^(r)^(r), (11) 

«W(r,t) = 5>£>(t)&(r)fc(r). (12) 



Real space non-condensate density and non-condensate correlations are, in general, nonlocal functions of two spatial 
points p(r',r) and «(r',r). We only computed these quantities for r = r' in this paper since these are the most 
physically accessible. Measuring these quantities with r^r' involves observing atomic correlations which is much more 
difficult than photon correlations. In Liouville space notation, the position-dependent nth order solution ip^- n \r,t) 
can be defined using the relations Eq. I|1(JI - I12|) and introducing a 2N(2N + 1) x 27V (27V + 1) square matrix T(r): 

</> } (r,£) = T(r)</»(t) (13) 

where 

T(r)=diag[0(r) ) ^(r),^(r),$ K (r),$;(r),$*(r)] . (14) 

Here "diag[- • •]" denotes that T(r) is a block diagonal square matrix made of TV x N blocks 4>(r), <p*(v) and N 2 x N 2 
blocks 4> p (r), <I> K (r), <j>*(r), <j>* (r). </>(r) is a diagonal matrix with the z-th diagonal element given by the basis states 
<j>i(r), and 3> P (r) and $ K (r) are also diagonal matrices whose zj'th diagonal element are given by [$ p (r)]^ ^. = 

4>*(r)4>j(r), and [$«(r)]y^ = 4>i{ v )4>] ( r ) respectively. The real space variables z^(r,t), p( n \r,t), and K^ n '(r,t) are 
finally obtained by summing over the appropriate elements of the vector ^iW(r,t): 

n 2n+n 2 2n+2n 2 

2W(r,f)=^f»(r,f), P (n) (r,t)= £ K (n) (r,i) = £ $ n) (r,t)- (15) 

i=l i=2n+l i=2n+ri 2 +l 

(2) 

The position-dependent second order (n = 2) response function K a (t, t±, t 2 , r, ri, r2) is then defined as follows: 

a (2) (r,<) = J K^[t,t 1 ,t2,r,r 1 ,r 2 ),Vf{r 1 ,t 1 )Vf(r 2 ,t 2 )dt l dt 2 dr 1 dr 2 , (16) 
where a^(r,t) are second order solutions, a = z,p,n. Expression for is given in Appendix IU1 Eqs. (|CJ1|) - (|CJ3|) . 

(2) 

Having found the time domain response K a (t, ti,t 2 ,r,ri,r 2 ), the second order susceptibility is obtained by a Fourier 
transform to the frequency domain: 

/>oo pOO />OC 

X^ ) {n,n 1 ,n 2 ,r,r 1 ,r 2 )= dt dh dt 2 K^{t,h,t 2 ,r,r 1 ,r 2 )eK.p(im + in 1 t 1 + iVL 2 t 2 ). (17) 

J o Jo Jo 

A closed expression for used in our numerical calculations is given in Appendix |Dl Eqs. l|DT) - (|D3l . 

III. TIME DOMAIN RESPONSE 

So far, all our results were given in the trap basis, and hold for a general interatomic interaction potential. In the 
following numerical calculations, we approximate the interatomic potential V(r — r') in Eq. JSJ by a contact potential: 

V(r-r')^U S(r-r'), U = ^^, (18) 

TO 

where a is the s-wave scattering length and to is the atomic mass. This is valid because wave functions at ultracold 
temperatures have very long wavelengths compared to the range of interatomic potential implying that details of the 
interatomic potential become unimportant. The tetradic matrices Viju defined in Eq. (JBJ are then simply given by: 

V ijU = / (r)^(r)^(r)&(r)dr. (19) 
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We consider a 2000 atom one dimensional condensate in a harmonic trap. The parameters used for our numerical 
calculation of ip^ are: U — 4ti ^' a = 0.01, and temperatures 0tki} tTa , p /k and 107iu; tra p/fc where w tr ap is the trap 
frequency, k is the Boltzmann constant and the basis set size of N = 5 is used. We keep the trap units throughout 
with 256 grid points for position. The same parameters were used in the calculations of the linear response in Ref. |l(. 

To solve for the second order response, both the zero'th and the first order solutions must be found. Calculation 
of the zero'th order solution from the TIHFB equations is the most numerically involved step, as it requires solving 
nonlinear coupled equations. Griffin has provided a self-consistent prescription for solving the TIHFB, in terms of 
the Bogoliubov-de Gennes Equations [6j . We have therefore followed the prescription of Ref. [6j to find the solution 
to TIHFB. Once the zero'th order solution is found, it is straightforward to calculate the first and the second order 
response functions. The calculation of the eigenvalues of the non-Hermitian matrix required for computing the 
response functions was carried out using the Arnoldi algorithm |3l| . 

In order to provide an indication of the structure of the matrix , we first plot in Fig. 1 the linear susceptibility 
JfW(fi,r = 0,ri = 0) for zero and finite temperatures. Peak positions indicate the resonant frequencies. 

An expression for the second order time domain response function is given in Appendix IO Eqs. I|C1|) - l|C3(l . 
We have first obtained the numerical solution to TIHFB, 2N(2N + 1) x 1 vector ipW evaluated at zero and finite 
temperatures, the 2N(2N + 1) x 2N(2N + 1) matrices f , U, and $ defined in Eqs. O, l|H6|l. and jOTfr. and the 
2N(2N + 1) x 1 vector E K defined in Eqs. (|C10IC13|) . Substituting these into Eqs. 1C4IC5|) . the final calculation 
involves matrix multiplication of 2N(2N + 1) x 1 vectors ^ (0) and E K with 2N(2N + 1) x 2N(2N + 1) matrices T, 
U, and and integration over the time variable r. T and $ are constructed in terms of the harmonic oscillator basis 
states which are calculated numerically from the recursive formula that involves the Gaussian function multiplying the 
Hermite polynomials [12 . The matrix U is calculated using a MATLAB function that uses the Pade approximation 
for matrix exponentiation |33| . 

We present the second-order response function in the time domain (t, t±,t2, r, ri^) as a function of r and ri at 
various times t, t±, ti and Vi- This provides a way to depict graphically the correlation involving 6 variables t, t\, £2, r, 
i"i,and Y2 on a 2 dimensional plot, and gives a "snapshot" of the position-dependent second order correlations across 
the condensate. The times i, t\, £2, are respectively the time of detection and the time of the first and second applied 
short fields, while r, ri,and r2 denote the corresponding spatial variables. The position-dependence is important 
since the experimentally produced condensates are mesoscopic in size; in optical spectroscopy, however, the dipole 
approximation usually applies and consequently the spatial dependence of the response is irrelevant. Fig. 2 shows 
the absolute value of the second order response functions in the time domain with the time ti fixed at £2 = 0. Fig. 2 
(a) is for the position of perturbation fixed at = i.e. the center of the trapped atomic cloud while Fig. 2 (b) is 
for T2 = —5 at the edge of the cloud. All positions are referred to in harmonic oscillator length units. The plots are 
for zero temperature condensate at the short, intermediate and long times ({t,ti} = {5.89,2.6}, {t,t{\ = {15.7,7.2}, 
{t, ti} = {31.4, 15.7}) as indicated at the top of each column of figures. The times are given in units of the trap period, 
l/wtrap- The top, middle, and bottom rows give the response function for the condensate z, non-condensate density 
/?, and non-condensate correlation k respectively. The dashed circle represents the spatial extent of the trapped BEC. 
The corresponding plot for a finite temperature BEC at temperature T — lOTnut^p/k is displayed in Fig. 3. 

The spatially asymmetric functions at r2 = — 5 reverses its shape for r2 = 5, with similar changes also observed for 
the Y2 = ±2.5 pair. With Y2 = ±2.5 which is a point inside the atomic cloud, the contours had markedly less symmetric 
shape than for Y2 — ±5 at the very edge of the atomic cloud. The response functions take a more symmetric shape 
as time increases. The Y2 dependence of the response function is found to be reduced at longer times. Comparing 
Figs. 2 and 3, the response functions clearly show strong temperature dependence with the functions giving distinct 
contours at different temperatures. The response functions attain spatial symmetry more rapidly at zero temperature 
owing to the weaker coupling between the variables z, p and k. 



IV. FREQUENCY DOMAIN RESPONSE 

Using Eq. (|DT4)l we have calculated the second order susceptibility. Eq. (|D14|I is also a matrix multiplication 
involving 2A(2A+1) x 1 vectors i/> (0) and E K (u) defined in Appendix|Dl Eq. (|D10)l . with the 2A(2A+1) x 2A(2A+1) 
matrices T, U(u>), and The Green's function U(uj) is calculated as follows: 

%) = — 1 L— = E 6< V ' ( 2 °) 

lu - £vA) -f te ^ u — u u + ie 

v 

where is the right eigenvector of with eigenvalues uj v such that C^^ y = lo v £, v and Cv are the left eig envectors 
of £^ 2 ) such that ^ y £„£t = 1. The eigenvalues uj v of were calculated using the Arnoldi algorithm 31] . 
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The absolute value of the second order response function in the frequency domain, \K^- 2 \fl, fix, fi 2 , r, rir2)| is 
displayed In Fig. 4 with the position variable r 2 set at (a) r 2 = at the center of the atomic cloud, and (b) r 2 = —5 
at the edge of the atomic cloud. The plots are for zero temperature condensate at various frequencies f^i and il 2 
indicated at the top of each column. We chose the frequencies such that f2i, £! 2 , and fix + f2 2 are off-resonant with 
respect to the eigenvalues of (first column, fix = 2.23, il 2 = 1.55); both fix and f2 2 are on-resonance while 
tlx + fl2 is off-resonant (second column, Qi = 2.2, £1 2 = 1.5); and finally Sli + f2 2 and 2 on-resonance (fix = 0.7, 
f2 2 = 1.5). As with the previous figures, the top, middle, and bottom rows give the response function for the 
condensate, non-condensate density, and non-condensate correlation respectively. The result for a finite temperature 
BEC at temperature T = lOTnot Yap / k is displayed in Fig. 5. At finite temperature, the resonant frequencies are 
shifted from the zero temperature counterpart so that the actual frequency combinations used are the following: Qx, 
f2 2 ) and + fl 2 off-resonant (Qj = 2.45, f2 2 = 1.6); both Six and f2 2 are on-resonance while fix + fl 2 is off-resonant 
(Hi = 2.43, S! 2 = 1.5); and frequencies with Sli + fl 2 is on-resonance (fix — 0.92, fl 2 = 1.5). 

The susceptibilities become more spatially symmetric as resonant frequencies are matched. The most symmetric 
function was generated when the sum of two frequencies fix and fl 2 was on resonance while the least symmetric 
function resulted when both frequencies were off-resonant. With only of one of the frequencies on resonance, a result 
with an intermediate level of symmetry was observed. The dependence of the response function on r 2 is strongest for 
the off-resonant case, while the one with fix +fl 2 on resonance remains more or less unaffected by the changes. As in 
the time domain response, the asymmetric functions at r 2 = —5 reverses its shape for r 2 = 5. Such change was also 
observed for the r 2 = ±2.5 pair; again, for r 2 = ±2.5 which is inside the atomic cloud the contours were found to be 
less symmetric than the corresponding figures for r 2 = ±5. 

In Fig. 6 we plot the second order susceptibility at zero temperature as a function of fix and fl 2 with the position 
variables set at r = ri = r 2 = i.e. at the center of the atomic cloud. The left column shows \K^ 2 \flx, ri 2 )| i.e. 
the absolute value of the second order response function. The plot shows only few peaks near the frequency of 1 
because the difference between the magnitude of these highest few peaks and the rest of the peaks occurring at other 
frequencies is too large. This suggests that one should ideally tune into the combination of frequencies fix and £1 2 
corresponding to these highest peaks to observe the maximum second order response experimentally. The second 
frequency fl 2 implies a completely different physics compared to the linear response; fl 2 is the frequency of a new 
harmonic being generated as a result of a strong external perturbation oscillating at frequency Q\. The middle column 
of Fig. 6 gives \K^ (fix, f2 2 )| where the largest peaks are scaled down to the magnitude of the smaller peaks present. 
It clearly shows a number of peaks present at frequencies f^i and fl 2 less than 1, and also shows that with fix (^2) 
fixed at 1, there is a pronounced response for many values of fl 2 (fix)- The right column of Fig. 6 shows the logarithm 
of the left column. This enables the large variations in the magnitude of \K^ 2 \flx, fl2)\ to be displayed. The top, 
middle, and bottom rows give the response function for the condensate, non- condensate density, and non-condensate 
correlation respectively. Fig. 7 shows the corresponding plot at the finite temperature T — lOfkutmp/k. The main 
difference is that there are more peaks appearing around frequency 0<O»<l,i=l,2 than at zero temperature. 

Since the second order susceptibility K^^flx, fl 2 ) is by itself a product of several Green's functions U(uj) and linear 
susceptibilities K^(uj), uj = fix, ^2, and fix + ^2, several features of Figs. 6 and 7 are closely related to the linear 
susceptibility. From Fig. 1, it is clear that the linear susceptibility near frequency of 1 dominates the spectrum at 
both zero and finite temperatures. It is therefore expected that all second order response K^ 2 \flx, SI2) containing 
f^i = 1 or fl2 = 1 component in it will lend a much stronger contribution to the spectrum than those not at fix = 1 or 
fl 2 = 1. This explains the features around fix = 1 and fl 2 = 1. In addition, the linear susceptibility K^(uj) at finite 
temperature has more peaks around frequencies < lu < 1 than at zero temperature, as there are more resonances 
below the frequency of 1 at finite temperature. This leads to the second order susceptibility K^(flx, fl 2 ) displaying 
more peaks around frequency < fii < 1 and < fl% < 1, as illustrated in Fig. 7. 



V. NONLINEAR RESPONSE OF FRACTIONAL QUANTUM HALL SYSTEMS 

The Fractional Quantum Hall Effect (FQHE) [2(j for a two-dimensional (2D) electron gas in a strong, perpendicular, 
external magnetic field is observed through the quantization of the Hall dc conductivity <jh{v) as a function of the 
filling factor v = (2-Kn)/(muj c ), where n is the mean 2D density of the electrons, lu c = eB/mc is the cyclotron 
frequency, B is the magnetic field; e is the electron charge (we set c = 1 and % = 1). At v = 1/(2 A; + 1), where k is 
an integer k = 1, 2, . . . <jh{v) varies in discrete steps, and is given by auiy) = v(e 2 /2tt). 

It has been shown that for v = l/(2k + 1) the original 2D fermion problem can be mapped into a boson system 
coupled to a Chern-Simons gauge field added to the time-independent external magnetic field , w here the original 
fermion system and the boson system have the same charge density and Hall conductivity [2l|, |22, UM ■ The static Hall 
conductivity (xy component of the conductivity tensor) calculated for this boson system, an = v(e /2tt), coincides 
with the conductivity obtained by using Laughlin's ansatz for the many-electron wavefunction |24|. For an even 
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inverse filling number v = l/(2fc) the original fermion problem is mapped onto a composite fermion system coupled to 
a Chern-Simons gauge field [25|]. At v = 1/2 this gauge field eliminates the external effective vector potential for the 
effective composite fermion system pfil 1271 l2*sj . In this section we consider FQHE with the odd inverse denominator 
of the filling number v = I /(2k +1). Using this mapping we can apply our results based on a generalized coherent 
state (GCS) ansatz for the many particle wave function of a weakly interacting effective boson system 0, 13M l36| to 
compute the second-order response of FQHE. Below we establish the correspondence between the FQHE effective 
Bose Hamiltonian and the Hamiltonian Eqs. Q- (J5J. 

We start with the many-electron effective Hamiltonian H where the electron-electron interaction contains Coulomb 
repulsion VJr — r') = l/|r — r'|, external magnetic vector-potential A(r), and an external electrical scalar potential 

v f (r,t) 

H = Jdt J dr^t) ( HVr " 2 ; (A(F)))2 -M) Mr,t) 

+ IJdtJdr J dr'i>i(r,t)^i(r\t)V(r-r')rl> e (r\t)i> e (r,t) 

dt /dr$(r,f)^(r,t)&(r,t), (21) 



where f/4(r, f) and ip e (r,t) are Fermi creation and annihilation operators; m is the effective electron band mass and 
/j, is the chemical potential. 

We next introduce a quasiparticle creation operator ffi(r,t) which is related to i/>|(r,i) as[26j 



ft(r,i) = ^|(r,t)exp 



J dr' arg(r — r')n(r', t) 



(22) 



where arg(r — r') is the angle between the vector (r — r') and the direction of the Hall current (which is perpendicular 
to both the external magnetic field B = V r x A(r) and the external electric field E = — V r V/(r, t), which are 
perpendicular to each other); v" 1 is an odd integer; and n(r,t) is the charge density operator, which is the same for 
the actual Fermi and artificial Bose system 

n(r,t) = $(r,t)&(r,f) = ft (r,t)rj;(r,t). (23) 

The operators ft(r, t) and ${r,i) satisfy Bose commutation relations |2lll22^ 

^(r)^t(r') - ft(r')4>{r) = S(r - r'); 
^(r)^(r') - ^(r')^(r) = 0; 

ft(r)ft(r') - ft(r')ft(r) = 0. (24) 

Using the Bose quasiparticle operators Eq. I)22|l , the Hamiltonian Eq. (|21JI can be written in a form of the many-boson 
effective Hamiltonian H by adding the gauge Chern-Simons vector potential a(r) pHl23| 

+ - f dt [ dr [ dr'^t( r>t )^t( r ' ( t)v(r-r')^(r',t)^(r,t) 



dt / drft(r,t)Vf(r,t)i>(r,t). (25) 



In order for the static Hall conductivity of the boson Hamiltonian Eq. 1251) in Hartree-Fock-Bogoliubov approxima- 
tion |3?| to coincide with the conductivity obtained by using Laughlin's ansatz for the many-electron wavefunction |24j , 
the magnetic and gauge potentials in k space should satisfy 23] 



A a (k)+a a (k) = — e a ^{h k -n), (26) 
ev k 



where e a ^ is a unit antisymmetric tensor 
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n is the mean 2D electron density n = (mvu) c ) /(2tt), and the Fourier components of the charge density operator fik 
with momentum k is 



= (28) 



where at and are Fourier components of f/^(r) and V'(r). Using a basis set of single electron functions </>i(r), the 
field operators may be expanded in the form 

^t( r>t ) = ^#( r )at(t); 

i 

^(r s t) = ^^(r)a i (t), (29) 

i 

where a] and are Bose operators. 

We make the following GCS Hartree-Fock-Bogoliubov ansatz for the time-dependent many-boson wavefunction 0, 

m 

= exp fX)oi(t)aI + j ( 3 °) 

where |fio) is an arbitrary normalized reference state with (Ool^o) = 1 |37| . 

On the other hand, the operator set of GCS generators of the closed algebra [3{| |3(| in the exponent of the 
ansatz Eq. I|30|l creates an extended Heisenberg-Weyl algebra which may be obtained by a repeated application of the 
standard boson commutators. 

The Hamiltonian that describes the system of many-body interacting bosons is obtained and comes from the basis 
set expansion of the Hamiltonian Eq. 12ol) and is given by [23] 

H = y^gjj-otoj + '^V i j k ia\a]a k ai + t] ^ (t)aj%, (31) 

ij ijkl ij 

Note that the Hamiltonian Eq. (|31J) coincides with Eqs. Q- JHJ) provided we replace in Eqs. (0 H*j by ffy. In 
Eq. (|3 1|) the single-electron matrix element, given by 

Hli= j dr ,. (r ^t^<m±^mi.,y l{t) , (32) 

VijA,; is a Coulomb repulsion between two electrons (we put electron charge e = 1) 

Vijki = [ dridr 2 </)*( r i) < / ) K I ' 2 )i " r0fe (ri)^ (r 2 ), (33) 

J J |ri — r 2 | 

and the external electrical field E(t) field expanded in a basis set is 

E^{t)a\a 3 = [ dr<t>*(r)V f (r, 0^(r)a]a,, (34) 

ij J 

We can, therefore, apply the results obtained for the Hamiltonian Eqs. (JIJ- J5j to the FQHE by simply replacing 
in Eqs. © by H t] Eq. 

Since in a homogeneous system the momentum is conserved, we can use the plane wave basis, i.e. the cigenfunctions 
of a momentum p (4> p (r) = exp(— ipr); where S is the 2D volume of a system) by replacing the indices i by 

the momentum p. 

The second-order response is given by Eqs. (Dl)- (D3),(D24) in the plane wave basis, provided we use for the 
Hamiltonian Eqs.(l)- (2) in the Liouvillian Eq. (B6) with Vt rap = 0, and replace the inter-particle interaction Vijki by 



27re 2 47r 2 ?i 

P_P ' = e|p-p'| + m^ 2 |p-p'l 2 ' (35) 



8 



where the first term in the r.h.s. represents the 2D Fourier transform component of the Coulomb repulsion and the 
second term comes from substituting Eq. i|2(ill in Eq. (|31[1 . With these substitutions, the FQHE boson Hamiltonian 
Ea (THTl reduced to the classical form Eq. (19) [23] and is equivalent to the classical boson Hamiltonian (Eq. (40) 
in [33). The double excitations are given by the eigenvalues of the Liouvillian Eq. (B6). The anisotropy in the gauge 
term Eq. (|26|l in the Hamiltonian Eq. (|31|l gives a non- vanishing second-order response, shifting the double-excitation 
energies with respect to twice the single excitations, as is the case in isotropic systems. 

VI. CONCLUSIONS 

We have calculated the second order mechanical response functions and susceptibilities of both zero and finite 
temperature BEC. The systematic, perturbative solution of the TDHFB equations for trapped, atomic BEC enables 
us to analyze the dynamics of finite temperature BEC. The response of both the non-condensate atoms as well as the 
condensate were calculated. The calculations apply for a general perturbation of the form J^ij I rf 3 r 0* (r)V/(r, t)a|dj , 
where the shape of the external force Vf(r,t) is left arbitrary. The calculated second order response functions were 
found to show a strong dependence on position and temperature. In the frequency domain, we have observed distinct 
nonlinear mixing effect, which displays enhanced response at a number of combinations. The application of 

the GCS ansatz for the derivation of the second-order response of FQHE 2D electron liquid was presented. 

The second order response for BEC has not been reported yet, and most related work is limited to the linear 
response. As in nonlinear optics, the second and higher order response functions and susceptibilities are clearly 
expected to be indispensable in characterizing the BEC dynamics in the presence of strong perturbations, and more 
importantly, in further development of matter- wave nonlinear optics 2] . 
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APPENDIX A: THE TDHFB EQUATIONS 

The TDHFB equations of motion that couples z, p and n are 0,0] 

dz 

ih— = [H z +r)E{t)]z + H z *z*, (Al) 

ih^ = [h, p] — (/tA* — Ak*) + r][E{t), p], (A2) 

ih^l = (hn + nh*) + ( P A + Ap*) + A + r)[E(t),K] + , (A3) 
dt 

where [. . .]+ denotes the anticommutator. Here, H. z , 7i z *, h, and A are N x N matrices with n being the number of 
basis wave functions used: 

= H H ~ M + Vikl * + 2pik ^ ( A4 ) 

kl 

[H z *] itj = yVijuKu (A5) 

kl 

hij = Hi-j - p + 2 ^2 Vikij \z* k zi + pik] , (A6) 

kl 

Ay = ^ V Vkl \ Z kZl + K k l] ■ (A7) 
kl 

h is known as the "Hartree-Fock Hamiltonian" and A as the "pairing field" ■ p is the chemical potential introduced 
in the Hamiltonian, Eq. J2J). 
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APPENDIX B: MATRIX C AND VECTOR C 



In this Appendix we define the 2A(2A + 1) x 2N(2N + 1) Liouvillian matrix of Eq. ©. As discussed in the main 
text, it suffices to define the zero'th and the first order matrices £C°) and only, since for all n > 1 is identical. 
First of all, £C°) is defined as follows: 



£ (0) = 



' riz — M 
















\ 


/Tz* 


Atz — /i 






















H(-) 







V 










£>a* 


WW 


V 
















p* 


n { -> 






V o 





D* 





_ V A 


w (+) * 


/ 



(Bl) 



where W z and TL Z * of are as defined in Eqs. I|A4IA5|) and the remaining A 2 x A 2 submatrices are defined as: 



ij,mn 




h (0) 6- 


ij : mn 


= fen + 






= a (0) <5 7 „ 






= -A* (0) <5 mi 



with and A(°) being the Hartree-Fock Hamiltonian and the pairing field defined in Eqs. 1A6IA7|) . 
For higher orders n > 1 (in particular, n = 2), 

where £( ) is the zero'th order matrix defined above, and 

)zzl ~\)zz2 ~\)zl ~\]z2 



(B2) 
(B3) 
(B4) 
(B5) 



(B6) 



a = 



I v zzL V 222 V 2i V 

yzz2* yzzl* g q 

ypzl ypz2 yyp'i W kA 

yKzi ynz2 yy^' 1 yy^^ 

ypz2* ypzl* g (yy/iAf-)* (yyph)* (yytcAj* 

V v Kz2 * y^zi* (yy«'it)* q (>V k ' 1 )* (W pA )* / 



\ 

W KAt 



(B7) 



The set of A x A submatrices V zzl and V 222 , A x A 2 submatrices V zl and V 22 , N 2 x N submatrices V p2 \ V p22 , 
V Kzl , and V K22 , and A 2 x A 2 component submatrices W ph , W pA , W Kh , and W kA of £' are given as follows: 



v zzi _ sr^ v *(o) jo) v zz2_sr^ v Jo) jo) 

V i,l ~ VtklrZ k Z r V i,k — 2-^i V *klrZ l Z r 



kr 



• z {0) 



(0) 



Zr lr 

v kz1 _ oV^x/ *(0) (0) T/ *(0) (0) [ T/ (0) T/ (0)' 

■^J.fc _ 2 / J VilkrZl K r j +VrkljZ l K ir + [ V rjklZ l + V r jlkZ l 

lr lr 



K 



K 



.(0) 
(0) 



+ > \Vtrklz\ + VtrlkZ\ 



*(0) 



y K z2 _ j^u 7 (0) (0) y (0) (0) 
~ Z / , v iklr z i K r j "I" Vrlk]Z l K, 



pA _ 



ij,kl 



^ K>fc;p r " * + V r jkip 



(0) 



(B8) 
(B9) 
(BIO) 
(Bll) 

(B12) 
(B13) 
(B14) 
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><fcj 



(0) 



r 



(0) 



reAf 



(0) 



(B15) 
(B16) 



In addition, the 2N(2N + 1) x 1 vector £(t) is defined as follows: 



C(*) = 






e(-)(t) 






V o 








(<) 











p<0) 



,(+) 



(*)]*/ 



g(o)« 



where -E(i) is as defined in Eq. JIJ, and we have further defined the iV 2 x N 2 submatrices 



(B17) 



(B18) 



APPENDIX C: SECOND ORDER RESPONSE FUNCTIONS 

In this Appendix, we first summaries the final result for the response functions, and then provide a more detailed 
derivation. In the last subsection of this Appendix, we also show an alternative form for the response functions written 
in the basis of the eigenvectors of the Liouvillian. 



1. Final expression 



The second order response functions for the condensate, the non-condensate density and correlation that we use in 
our numerical calculations are given by Eqs. IjClfl - (|C3I) as follows: 



N 

(t, t lt t2, r, n, r 2 ) = J2 [T(r) (itf \t, t x , ta, r x , r 2 ) + Rf} (t, t lt t 2 , r x , r 2 ] 



(CI) 



/fW(t,ti,t2,r,ri,r 2 ) 



2N+W 



i=2N+l 



[f (r) (Rf ] (t, h , ta, r x , r 2 ) + (t, ti , ta, r x , ra) 



(C2) 



2N+2N' 

lfW(t,tl,t2,r,ri,r 2 ) = [*( r ) (^f Wi,*2, r x , r 2 ) + X^V, i 2 , r x , r 2 ) 

i=2JV+JV 2 + l 



(C3) 



where 2N(2N + 1) x 2N(2N + 1) matrix f (r) was defined in Eq. (JUJ, and the 2N(2N + 1) x 1 vectors i?f ) and 

(2) 



Kjj are defined as follows 



^ 2) (t,ti,t 2 ,ri,r 2 ) =W(t-ti)$(r 1 )W(ti-t 2 )#(r 2 )V 



•(0) 



(C4) 



^}j ) (*,ti,t 2 ,ri,r 2 ) = / drW(i-r)Hx(r-ti,r-*2,rx,r 2 ), 



(C5) 



where 



W(i - f') = exp 



-C (2 \t-t') 



(C6) 
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and 



$(r') = diag f$(r'),$*(r'),$(-)(r') 1 $W(r / ),$(-)*(r / ),*( + >(r / ) 



(C7) 



$(r') is therefore a 2N(2N + 1) x 2N(2N + 1) block diagonal square matrix with the blocks consisting of N x N 
square matrices 



and N 2 x N 2 square matrices 



$(±)( r ) 



[$(r)].,=0*(r)^-(r) 



</>*(r)</> m (r)^„ ± ^(r)^-(r)4 



The vector H# of Eq. ljC5(l may be written 

with the TV x 1 matrix Zk and N 2 x 1 matrices TZk, ICk given as follows: 
[Z K \ = £V ifcIr [{Ki 1] *(T -t^r^z^lKPir -t 2 ,r 2 )] r 

klr 

+ [KP*(T - h^MK^iT - f 2 ,r 2 )]4°) + z^WV - ^nMi^r - i 2 ,r 2 )] T 



(C8) 



(C9) 



(CIO) 



+ 2[Kf\r - ti,n)], fc [lirW(T - i 2 , r 2 )] r + [tf W(t - ti,n)] w [tf W(t - fc, r 2 )] r 



«/ 



Hi)/ 



(Cll) 



rfe/ 

- V^y (r - t 1} r x )] fci (r - i 2 , r 2 )] ir + £ V**, [K^ (r - *i, n)] fcI (r - t 2 , r 2 )] f 



rA:/ 



+ V^[i?* (1) (V - ii.ri)]*,^ 1 )^ - i 2 ,r 2 )] lr 

+ 2^T/ 4feir [i?;W(r-i 1 ,r 1 )] fc [i?i 1) (r-t 2 ,r 2 )] i 4 , ; ) 



rA'/ 



2$> iHr [[^;^(r - ti )ri )] fc z< 0) + ^ ( °W } (t - ti,n)],j [i?«(r ~ t 2l r 2 ) 



rkl 



Vrjki [[2?P(r - ti,n)] k <W + z; (0) K«(r - fl) n)] z ] [i?W(r - i 2 ,r 2 )], 

V irU [[K?\t - tx.n)]*^ + - tx.ri)],] [#:«(r - * 2 ,r 2 )] rj 

X)^fcj[^ 1) (r-ti 1 r 1 )] fc [^W(r-t2,r 2 )] l «^ 



rA'/ 



<0) 



(C12) 



— t%, Yi)]ki[K^ (r — t 2 ,r 2 )]i r + Virfc/ [K^\t — t\, T±)]ki [KjP*(r — t 2 , r 2 )] r 



rA:/ 



- il,ri)] w [^^(r-t 2 ,r 2 )] ir + 2 ^ ^ fcir [7?* (1) (r - i l5 ri)] fe [i?W(T - i 2 ,r 2 )];rc 



rA.7 



+ 2^VW [[^ 1 )*(T-ti,r 1 )] fc z J (0) +«*W[^ 1 )(r-t 1) r 1 )] I ] [K^{r - h, r 2 )] r 



rkl 
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+ V rklj [[^(r-iLrO^'+^f [^(r-ti.rx)],] [i?«(r - i 2 , r 2 )]„ 

+ J2 V rM [[^ 1) (^-*i ) r 1 )]^ (0) + 4° ) [^i 1) (T-*i,ri)] i ] [xW(r-t 2 ,r 2 ); 

+ VS rW [[^(r-ii.rOl^ + ^^Cr-ti.r!)],] [>f;«(r - * 2 , r 2 )] rj 
+ V nkl [KP (r-h, rrMKP (r-t 2 , r 2 )]^ r 0) 



rkl 



<0) 



+ VirM [KP(T - h, nMKP (T-t 2 , V2)]ip\ 



(C13) 



where the vectors Kz(t — ti,Ti), Kp(t — ti,rj) and Kk(t — ti,Ti) found in Eqs. (|C11IC13|) are defined as 

itM(T-U,Ti) =U^(r-U)^(r i )^°\ a = z,p,K, i = 1,2. (C14) 
£/0)( r _ ^) are the submatrices of U(t — t\) such that 

U{t - h) = \u {z \t - U),U (z) *{r - U),U (p) {r - U)M k) (t - U)M (p) *{t - U),U {k) *{t - U)\ T , (C15) 

where U {z \t - U) is an N x 2N(2N + 1) submatrix while U^ir - U), 7 = p,K,p*,K* is an N 2 x 2N(2N + 1) 
submatrix, i.e. the submatrix U^ z \t — ti) is stacked on top of submatrix U^*[r — ti) which, in turn, is stacked on 

top of submatrices U^\t — ti). It is to be noted that Ka\r — U, r^) as defined here are N x 1 and N 2 x 1 vectors, 
not scalar quantities obtained by integrating the scalar function K„\t — U,r, r,) over r. 



Derivation 



Writing the second order solution to TDHFB explicitly, we have 



^ 3 >(t) = 



k 1 (,xp 



r(ti)dti 



— / U{t-ti)T{ti)dt x 



where 



(C16) 
(C17) 

(C18) 



i.e. for the second order response, \W{t) = C(*)$ (1) (*) +H(t) in Eq. © where f(t) is given in Appendix FBI and H(t) 
is a 2N(2N + 1) x 1 vector originating from the terms in the expansion which are made up of products of two first 
order variables i.e. z (1) , and The vector S(t) can be written as 3(t) = [Z, Z*K, JC, K*,1C*} T with the N x 1 
matrix Z, and N 2 x 1 matrices 1Z, JC given as follows: 



Z. 
K 



klr 

rkl rkl 

2fe ' r fc i "rj Vrkl 3 Z k Z l Pir 



(C19) 



r/,7 



2 ^ Vifeir 
rkl 



rkl 



*(1) JO) , *(0) (1) 



*(1) *(0) . *(0) *(1) 



P r .j - Vrklj 



(1) 



K 



VirM 



,*(1)„(0) , *(0) (1) 



-(l)-(0) , JO) Jl) 
^fe "T" z k z l 



S 1 ) 



,(1) 



13 



Y, V rjkl4 



(i) Ji)jo) 



Virkl Zfa 



(1) Ji) *(0) 



(C20) 



rkl 



(1) Jl) 



*(i)„ (1)^(0) , \t Ji) *(i) to) 



rfcZ 



*(1) JO) , *(0) (1) 
z fc *l z fc z i 



Jl) *(0) (0) *(1) 

z fc z fc z i 



.(1) 



rkl 



;A7 



Jl) JO) , JO) Jl) 
z fc z i z fc z i 



Jl) JO) , JO) Jl) 
z fe z / ' z fe z i 



IT] 



v^t/ _(i)_(i;jo T/ (l) i *(0) r: 
2_^V rj kiz k z / P lr +V irk iz k z l p rj + 2_^ v W z k 



Ji) Ji) 



rkl 



kl 



Casting Eq. (|C17|I in the form 

J (2) (r,t) = J K {2 \t 1 t l ,t 2 ,r,r l v 2 )V s {v ll t l )V } {r 2 M)d i Y 1 dt l d i Y 2 dt 2 
involves rewriting T(t) of Eq. (|C18|I in the position dependent form: 

$(n)i?« (h ,t 2 , r 2 )] Vf (ri, ti)Vf (r 2 , t 2 ) 

- *2,*i ~ £3;r 2 ,r 3 ) V/(r 2 ,t 2 )V/ (r 3 ,t 3 ), 



r(t x ) = ^ J *a / rfr i rfr 2 



o Jo 



dt 2 dts / dr 2 dr. 



where i?^ 1 -* (ti, t 2 , r 2 ) is the linear response function for the combined variables z, p and k: 

K {1) (ti - ta, r 2 ) = W(ti - i 2 )$(r 2 )^ (0) , 



(C21) 

(C22) 

(C23) 
(C24) 

(C25) 



and the 2N(2N + 1) x 1 column vector Sif(ti — f 2 ,ti — i3;r 2 ,r 2 ) has already been defined above in Eqs. IjClOt - 
IC13II . S# is derived from the vector S(i) such that the components of the linear response vector iCW(ti ,t 2 , r 2 ) i.e. 
J^i 1 - 1 (ti , i 2 , r 2 ), ^p 1 ''(ti,t 2 ,r 2 ) and K^\t\, t 2l r 2 ) defined in Eq. l|C14p replace, respectively, z^\ti), p^(t±), and 
KW(ti) in S(t). 

Using these results, the position-dependent time domain second-order response function for the combined variables 
z, p, and k may finally be written as 



K( 2 )(t,t 1 ,t 2 ,r,r 1 ,r 2 ) = f(r) 



where, after the change of variables t\ — » r, f 2 — > ti, t% — > i 2 , A") and if}/ are given by 

^i 2) (t,ti,t2,ri ) r 2 ) = W(t - ti)$(ri)W(ti - t 2 )$(r 2 )^°>, 



^j ) (t,ti > t 2> ri I r 2 ) 



drW(i - r)3 K (T -fi,T- t 2 ,ri,r 2 ). 



(C26) 

(C27) 
(C28) 



The position-dependent second order response functions for the condensate, the non-condensate density and corre- 
lation are given by summing over appropriate indices of the vector K^ 2 \ Eq. (|C26|I : 



(C29) 



N 

Kf> (t, t ! , t 2 , r , n , r 2 ) = Y, ( r ) (Ki 2 } (* , *i , *a , r x , r 2 ) + (i, t x , t 2 , n , r 2 ) 



^(t.ti.ta.r.ri.ra) = 



2N+N 

E 

i=2JV+l 



[f (r) ) (t, ti, t 2 , r x , r 2 ) + K.f] (t, t x , t 2 , n, r 2 ) 



(C30) 
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2N+2N Z 



KP (t, t u t 2) r, n, r 2 ) = (^i 2) (*» *i> *2, ri, r 2 ) + xg 3 (i, fa , i 2 , n, r 2 ) 



i=2iV+AT 2 + l 



(C31) 



wher e iff } an d are defined in Eqs. (TU27|) - IU28|) . and f , l> and S K are denned in Eqs. O, JHSj), (fUT^ . 
and HC10IC13P respectively. 



3. Alternative form for the response function 



To discuss the response functions in the frequency domain and to understand the physical processes involved, it is 
useful to expand the response functions in the basis of the eigenvectors £„, of matrix such that £( 2 '£„ = <j} v £, v , 
z/ = 1,2, . . . 2N(2N + 1). We define the Green's function 

GJt-t') = exp ~wM-i!) , (C32) 
n 

and the expansion coefOficients r) u (r), and S u (r) such that 

2AT( 2 Af+l) 2JV(2jV+l) 2AT(2Af+l) 

^ (0) = E k*)&= E i»(. r )£» and t»£ = E *"( r )^- ( C33 ) 

z^— 1 z^— 1 u— 1 

In the basis of these eigenstates, £„, the result is Eqs. (jC29() - ijC31|) but with 

2^(2^+1) 

Kf\t,t x ,t 2 ^v 2 ) = Y, ^S(*,*i,*2,ri,r a )^, (C34) 

2^(2^+1) 

^ ) (t,ti,*2,r 1 ,r 2 ) = ^ 4 2 i(^*l.*2,ri,r 2 )|;, (C35) 



where 



and 



2JV(27V+1) 



/C} 23 (<,fa,i 2 ,ri,r 2 ) = 53 ^( r i)^'( r 2)G,y(t - ii)/i„"Gv'(ti - t 2 



j/',i/"=i 



2JV(2JV+1) t 



^X/i/(*i*ij*2,ri,r 2 ) = 51 / drG u{t - T)T[f] l/l {v l )-q vll {Y 2 )^^ vl >G^ 
v>y=i Jo 



(r-fa)G^(r-t 2 )] 



(C36) 



(C37) 



where J 7 is the function given according to the expression for S# but written in terms of \i v , r\ Vl and G v with G„(i— £'), 
\i V i and ?7„ as defined in Eqs. l|C32IC33j) . 

APPENDIX D: SECOND ORDER SUSCEPTIBILITIES 
1. Final expression 



The final result that we use for our numerical calculation for the condensate, non-condensate density and the 
non-condensate correlations are: 

N 

tfWt-fii-fiajni.fia.r.rx.ra) - Y [*W ^2, ri, r 2 ) + ^(^i, fi 2 , r x , r 2 ))] ^ (Dl) 



2N+N Z 



Kf\-^ - fi 2 ; fix, fi 2 , r, r 1; r 2 ) = 53 [f (r) >(«!, fi a , n, r 2 ) + ^(fti, fia, n, r 2 



i=2N+l 



(D2) 



2N+2N* 



Jr( 2 >(-fii-n2;fii s n2,r,ri,r 2 ) = E [f(r)(^f ) (ni,fi 2! r 1 ,r 2 ) + ^ ) (n 1 ,Q 2 ,r 1 ,r 2 )^ 



=2N+N 2 + 1 



where 



and 



1 



iff '(fli, 2) r, n, r 2 ) = -_T(r)W(ni + fi 2 )$(r 1 )W(fi 2 )l>(r 2 )v7 (0) , 



8ttH 



T(r)W(n 1 +n3)Hjf(ni,f23,ri,r 2 ). 



Here, T(r) and 3>(r) are as denned in Eqs. I|14|l and (|C7|I and 



1 



w-£( 2 )+ie' 

In addition, the vector Sj<-(f2i, f2 2 , ri, r 2 ) of Eq. I|D5(I is may be written as 
with the N x 1 matrix Zk and TV 2 x 1 matrices TZk, ICk given as follows: 



Z K 



K 



K 



$>Hr [[^i 1) *(ni,r 1 )] fc * I (0) [^ 1 )(n 2) r 2 )] r + [K£>{n 1 ,v 1 )] k [K£\n 2 ,v 2 )} lZ W 

klr 

2 E [#< 1} (fii , rr)] ri (0 2 , r 2 )] 

rhl 

E ^ ("i. r i)W r 2 )] rj + v rjkl (Oi, n)]« [#< 1} (n 2 , r 2 )] ir 

rfci 

2XV <Hr [ir;< 1 >(fii,r 1 )] / ^ 
rfei 

2 E ^kJr r i)]fe^ 0) + ^ (0) («i> ( fi 2, r 2 )] rj 

V pWj [[^(nx.n)]^ +z* (0) [^ 1 )(fi 1) r 1 )] I ] [/?«(n 2 ,r 2 )]ir 
E^-i [[iff'*(!l ll r 1 | Z ; (0) + z; (0 »K( 1 )(!! ll r 1 )] 1 ] (fi 2 , r 2 )] ir 

rfci 

Vfrw [[^(ni.n)]*^ +4 0) [^W(n 1 ,r 1 )] z ] [^*W(f2 2 , r 2 )] r j 

E^i^^ri)]*^ 1 ^ 
rfci 

2 E Vikir (Oi, n)]rj (fia, r 2 )]j fc + V rWj (Six, ri)]w [i?^ (0 2) r 2 )] ir 

rfei 

E ^frfcl [^H^i. i-OMiTpC^, r2 )] rj + y rifei [^W(ni,ri)] H [^W(n2,ra)] ir 

rfei 

2EW^ (1) (^i,ri)M^ 1) (^2,r 2 ^ 

rfei 

2E^ [[^(ni.ti)]^^ + zf ) [KP(Sl 1 ,r l )} l ] [K^iSlx,^)}, 



rkl 



V rklj [[KPiSlr^)}^ +^\K*W(Sl 1 ,r 1 )} l ] [KP(Sl 2 , r 2 )], 
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Y,Vrjkl [[K^i^rt)]^ +zf ) [K| 1 )(n 1 ,r 1 )],] [J?«(Q 2 ,r 2 )] ir 

rkl 

V irkl [[Ki 1 \Q 1 ,v 1 )] k zl 0) +4°>[^W(n 1 ,r 1 )] l ] [^;W(n3,r a )] Pj 

X] Kifci[i?f 1) (^i,ri)]fc[i?f 1) (^2,i"2)]iP-° ) + Vi r fe;[i?f 1) (fii,ri)] fc [i?f 1) (fi2,r 2 )]ip^ 0) 

^ [^S 1 ) (O x , ri )] fe [^l 1 ) (0 2l r 2 )] 2 . (D10) 



kl 



Similar to Eq. I|C14|I above, the quantities Kl J (f2 f ,ri), ^ r^) and /ft '(Qi,ri) used in Eqs. HD8ID10|) are 
denned as 



where a = z,p,K, and i = 1,2. W( a )(f2j) are the submatrices of U{£li) defined in Eq. 1D6|I such that 

u(Qi) = \u {z \n t )M (z) \^)M (p \^)M (K \^)M (p) *mM (K) *m 



(Dll) 



(D12) 



where i/ (z) (Qj) is aniVx2iV(2iV+l) submatrix while W (7) (^), 7 = p,n,p*,n* is an TV 2 x 2iV(2iV+ 1) submatrix such 
that the submatrix (f2j) is stacked on top of submatrix W' z )*(f2i) which, in turn, is stacked on top of submatrices 
U^^i). ft is to be noted that, as with the time domain example discussed above, Ka\^h,^i) as defined here are 
N x 1 and jV 2 x 1 vectors, not scalar quantities obtained by integrating the scalar function Ka (fij, r, rj) over r. 



2. Derivation 

The second order response function in frequency is given by the Fourier Transform of the time domain counterpart: 



K^\n,n 1 ,n 2 ,r, n,r 2 ) 



dtdtiT(r) 



T(r) Kf \n, Q lt Q 2l ri,r 2 ) + K% J (n, flx,^2, n,r a 



exp (iftt + ifiiti + 15^2^2) 



(2), 



(D13) 
(D14) 



where Kj 2 \$ l, Q\, fl 2 , r 2 ) and Kj Z j(fl, fii, f2 2 , ri, r 2 ) are the Fourier transforms of the time domain expressions Eq. 

E23- i|c28|) . 

Using the fact that the matrices U(t — ti) are the Green's functions with an implicit Heaviside step function in time 
i.e. U(t - t x ) = 9(t - t x )U{t - tx) such that 



H2), 



9{t)U(t) 



1 

2tt7 



1 



w - £( 2 > + ?< 
dtuU(uj) exp(— iujt) 



— exp(— itot) 



(D15) 
(D16) 



we have 



Kf\n,nx,n 2 ,v,vx,v 2 ) 



1 



dwdw* f (r)W(w)$(ri)W(w')*(r2)^ 0) <y(n - w) 



x5(fti + w - w')<J(^2 + w'). 
This implies that u/ = — f2 2 , u; = —fix — ^2, and = — 17i — f2 2 : 



1 



'(0) 



4tt 2 



(D17) 



(D18) 



In addition we have, 



Kf}{n^x^2,v,vx,v 2 ) = - ' 



87T 3 i 



dojdoj'duj" f(T)U(u))§ K (u)\ uj", n, r 2 ) 



(D19) 
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We are able to write the Fourier Transform for Sjf(t) in Eq. (|D19|I since the function S#-(i) is made up of terms 
which are simply products of two Green's functions at different times. Eq. 1L) 19(1 implies u' = — fii, uj" = — S7 2 , 
uj = u>' + uj" = —Q,i — O2, O — lo so that 



KfX-O x - 2 ; O x ,0 2 , r, n, r 2 ) = --_T(r)W(fii + fi 2 )H K (fii, fi 2 , ri,r 2 ), 



(D20) 



where S^-(fii, fi 2 ,ri,r 2 ) is as already given in Eq. (|D10J| . 

As for the time domain calculations, the susceptibilities for the condensate, non-condensate density and the non- 
condensate correlations are obtained by summing over the appropriate indices: 

N 

K^(-0 1 -0 2 ;O u 2 ,v,r 1 ,r 2 ) = £ [f (r) (i?f \o 1: 2 , n, r 2 ) + ^(fii, 2 , r x , r 2 ))] . (D21) 



2N+N 



JfW(-fii - fi 2 ; fii, fi 2 , r, n, r 2 ) = £ [f (r) } (fii, fi 2 , n, r 2 ) + ffg^i. fi 2 , n, r 2 )) 



i=2N+l 



(D22) 



2N+2N 



K^(-0 1 -0 2 -0 1 ,0 2 ,v,v 1 ,v 2 )= [T(r)(i?f ) (r! 1 ,r! 2 ,r 1 ,r 2 )+^ ) (fi 1 ,r! 2 ,r 1 ,r 2 ))] i (D23) 

i=2Ar+Ar 2 + l 



3. Alternative form for the susceptibility 

As before, expanding in the eigenstate basis £„ which we introduced in Eq. l|C33p . we may write the susceptibility 
m a more useful form. The result is Eqs. l|D2Tj) - l|D23)l but with 



2Af(2Af+l) 



^ 2) (-n 1 -f2 a ;fii,n a ,ri,r 2 ) = £ - ^2; fii, 2 , n, r 2 )£,, 

2^(2^+1) 

^(-ni-fJajni.fia.ri.ra) = ^ - ^2; fii, fi 2 , r x , v 2 %, 



where 



2N(2N+1 

ICfl(-0 1 -0 2 ;0 1 ,0 2 ,v 1 ,v 2 ) = --^ 



,~f_ x (Oi + 2 - uj v + ie)(0 2 - uo v > + ie) 



(D24) 
(D25) 

(D26) 



(2) 

Since no additional information is gained by listing all the terms in Kjj v (—0\ — fi 2 ; fii, 2 , ri, r 2 ), we simply note 

(2) 

that, in the eigenstate basis, the typical term in /Cjt v has the structure: 



2N(2N+1) 



E 



^(ri)fy/(r 2 )/V< 



,„ = , (^1 + 2 - u v + ie)(0\ - u v + ie)(0 2 - u) v > + ie) 
as to be expected from Eqs. (|D10|) and (|D20(I . and rj v are as given in Eq.JU33j|. 



(D27) 
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Fig. 1: Natural log of linear susceptibility K^\0,r, ri) at r = ri = vs. frequency. Top three panels - zero 
temperature, bottom three panels - finite temperature 10huj/k. The frequency Ox are given in units of the trap 
frequency. 



Fig. 2: \KW(t, ti, t 2 , r, r! r 2 ) | i.e. the absolute value of the second order response functions in the time domain with 
the time t 2 fixed at t 2 = 0. (a) r 2 = 0; (b) r 2 = —5. The plots are for zero temperature condensate at the short, 
intermediate and long times t and t\ written at the top of each column of figures. The top, middle, and bottom rows 
give the response function for the condensate, non-condensate density, and non-condensate correlation respectively. 
The diameter of the dashed circle represent the spatial extent of the trapped BEC. The position x and x' are given 
in harmonic oscillator length units. 



Fig. 3: Same as in Fig. 2, but at finite temperature of lOfuv/k. The position x and x' are given in harmonic oscillator 
length units. 



Fig. 4: \K^ 2 \0, Q,i, 2 , r, i"ir 2 )| i.e. the absolute value of the second order response functions in the frequency 
domain with the variable r 2 set at (a) r 2 = and (b) r 2 = —5. The plots are for zero temperature condensate at the 
frequencies Oi and 2 given at the top of each column. Denoting "off-resonant" when the frequency dos not match 
an eigen value of and "on-resonance" when a frequency matches an eigen value, these frequencies are chosen 
such that fii, O2, and 0\ + 2 are off-resonant (0\ — 2.23, 2 = 1.55); both Oi and 2 are on-resonance while 

01 + 2 is off-resonant (f2i = 2.2, 2 — 1.5); and finally frequencies chosen so that f^i + 2 is on-resonance {0\ = 0.7, 

2 = 1.5). The top, middle, and bottom rows give the response function for the condensate, non-condensate density, 
and non-condensate correlation respectively. The diameter of the dashed circle represent the spatial extent of the 
trapped BEC. The position x and x' are given in harmonic oscillator length units. 



Fig. 5: Same as in Fig. 4, but at finite temperature of lOfkv/k. At finite temperature, the resonant frequencies 
are shifted from the zero temperature counterpart so that the actual frequency combinations used are the following: 
0\, 2 , and Oi + 2 off-resonant (Oi = 2.45, 2 = 1.6); both Oi and 2 on-resonance with 0\ + 2 off- resonant 
(0\ = 2.43, 2 = 1.5); and finally Oi + 2 on-resonance (0\ = 0.92, 2 = 1.5). The position x and x' are given in 
harmonic oscillator length units. 



Fig. 6: Left column: \K^ 2 '(Oi,0 2 )\ i.e. the absolute value of the second order response functions in the frequency 
domain as a function of Oi and 2 with the position variables set at r = ri = r 2 = 0. There are three large peaks 
that completely dominates the plot so that the remaining peaks are not represented in the plot. Center column: 
Scaled \R( 2 \Oi, 2 )\. The three largest peaks shown in the left column were scaled down to the same magnitude as 
other peaks in the plot. Right column: logji^^ 2 ^ 1 (Oi, f2 2 )| i.e. logarithm of the response functions presented in the left 
column to help visualize the large variation in the magnitude of \K^(Oi,0 2 )\. The top, middle, and bottom rows 
give the response function for the condensate, non-condensate density, and non-condensate correlation respectively. 
The frequency Oi and 2 are given in units of the trap frequency. 
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Fig. 7: Same as in Fig. 6, but at finite temperature of lQTwj/k. The frequency Oi and f2 2 are given in units of the 
trap frequency. 
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